#This to obtain:
#Figure_6b.jpg


rm(list = ls())
library(sf)
library(stargazer)
library(ggplot2)
library(dplyr)
library(lemon)
library(lfe)
library(broom)
library(tidyr)
library(lemon)
library(ggplot2)
library(grid)
library(reshape2)

#Step 1: Read Data
#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

census_1895_dist <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                            layer="census_1895_points")

census_1895_dist<-st_drop_geometry(census_1895_dist)
census_1895 <- subset(census_1895_dist, inside_croatia == 1)
names(census_1895_dist)
rm(census_1895_dist)

#Step 2: Perform calculation
#census_1895$treat<-census_1895$inside_krajna6_treat
census_1895$dist1 <- as.numeric( with (census_1895,ifelse(census_1895$treat==1, 1, -1)))
census_1895$krajna6_distance<-census_1895$krajna6_distance/1000
census_1895$zagreb_distance<-census_1895$zagreb_distance/1000
colnames(census_1895)
data <- subset(census_1895, treat == "1" | treat == "0")
data<-subset(data, area!="*)")

data$dist2 <-data$dist1*data$krajna6_distance
data$bfe1 <- ifelse(data$krajna6_NEAR_FID == 1, 1,0)
data$bfe2 <- ifelse(data$krajna6_NEAR_FID == 2, 1,0)
data$bfe3 <- ifelse(data$krajna6_NEAR_FID == 3, 1,0)
data$bfe4 <- ifelse(data$krajna6_NEAR_FID == 4, 1,0)
data$bfe5 <- ifelse(data$krajna6_NEAR_FID == 5, 1,0)
data$bfe6 <- ifelse(data$krajna6_NEAR_FID == 6, 1,0)
data$bfe7 <- ifelse(data$krajna6_NEAR_FID == 7, 1,0)
data$bfe8 <- ifelse(data$krajna6_NEAR_FID == 8, 1,0)
data$POINT_Y<-data$lat
data$POINT_X<-data$lon

data$quad1<-ifelse((data$lon > 18 & data$lon< 20) & 
                     (data$lat<47  & data$lat>44), 1, 0)


data$quad2<-ifelse((data$lon > 16 & data$lon< 18) & 
                     (data$lat<47  & data$lat>44), 1, 0)


data$quad3<-ifelse((data$lon > 14 & data$lon< 16) & 
                     (data$lat<47  & data$lat>44), 1, 0)


data$log_existing_zadrugas_yokes<-log(data$existing_zadrugas_yokes+1)


m1 <- felm(log_existing_zadrugas_yokes~treat+
             POINT_X+POINT_Y+zagreb_distance+
             bfe1+bfe2+bfe3+bfe4+
             bfe5+bfe6+bfe7+
             quad1+ quad2
           |0|0|county, 
           data = data)
summary(m1)
m1<-tidy(m1)




#This saves all the previous regression results in a dataframe.
all_models <- bind_rows(
  m1 %>% mutate(model = "log_existing_zadrugas_yokes"))

#This saves the OLS mode by gathering variables.
ols_table <- all_models %>%
  dplyr::select(-std.error, -statistic, -p.value) %>%
  #mutate_each(funs(round(., 2)), -term) %>% 
  gather(key, value, estimate)
#spread(model, value) 


#The variables become column names.
#Models become rows.
regcoef<-dcast(ols_table, model~term, fill=0)
colnames(regcoef)
colnames(regcoef)[2] <- "Intercept"


########
#Part 2#
########
#load packages
library(foreign)
library(ggplot2)
library(rgdal)
library(RColorBrewer) # creates nice color schemes
library(maptools)     # loads sp library too
library(scales)       # customize scales
library(gridExtra)    # mutiple plots
library(plyr)     
library(dplyr)
library(grid)
library(tidyr)
library(ggsn)
library(lemon)



#### Load data that will be used in all the figures

boundary <- st_read(dsn='./Dropbox/Legacies_Project/Analysis/data/data.gdb',
                    layer="krajna_line_HRV")


grid <- st_read(dsn='./Dropbox/Legacies_Project/Analysis/data/data.gdb',
                layer="grid4c")


boundary <- st_transform(boundary, crs = st_crs(grid))


#colnames(grid.df) <- c("Id", "treat", "x","y", "dist_bnd", "bnd", "dist_hcmc")
grid_df <- grid%>% mutate(zagreb_dist = (dist_zagreb/1000),
                             krajina6_dist = (dist_krajna6/1000))

grid_df$quad1<-ifelse((grid_df$xcoord > 18 & grid_df$xcoord<20) | 
                        (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)

grid_df$quad2<-ifelse((grid_df$xcoord > 16 & grid_df$xcoord<18) | 
                        (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)

grid_df$quad3<-ifelse((grid_df$xcoord > 14 & grid_df$xcoord<16) |
                        (grid_df$xcoord<47 & grid_df$xcoord>44), 1, 0)

unique(grid_df$quad3)

#The next line is to get border FE
grid_df$krajina6_id1<-0
grid_df$krajina6_id1[grid_df$krajina6_id==1]<-1
grid_df$krajina6_id2<-0
grid_df$krajina6_id2[grid_df$krajina6_id==2]<-1
grid_df$krajina6_id3<-0
grid_df$krajina6_id3[grid_df$krajina6_id==3]<-1
grid_df$krajina6_id4<-0
grid_df$krajina6_id4[grid_df$krajina6_id==4]<-1
grid_df$krajina6_id5<-0
grid_df$krajina6_id5[grid_df$krajina6_id==5]<-1
grid_df$krajina6_id6<-0
grid_df$krajina6_id6[grid_df$krajina6_id==6]<-1
grid_df$krajina6_id7<-0
grid_df$krajina6_id7[grid_df$krajina6_id==7]<-1
grid_df$krajina6_id8<-0
grid_df$krajina6_id8[grid_df$krajina6_id==8]<-1

regcoef[is.na(regcoef)] <- 0


#################
#log_existing_zadrugas_yokes#
#################

summary(data$log_existing_zadrugas_yokes)

subgrid <- grid_df
subdata <- data

subdata <- subset(subdata, krajna6_distance<=30)
subcoef <- subset(regcoef, model=="log_existing_zadrugas_yokes")
subgrid <- (subgrid  
            %>% mutate(ones=1))  # Create a column of ones for constant


subgrid <- data.frame(lapply(subgrid, function(x) as.numeric(as.character(x))))


subgrid$fit1<-subcoef$treat * subgrid$treat+
  subcoef$POINT_X * subgrid$xcoord+
  subcoef$POINT_Y * subgrid$ycoord+
  subcoef$zagreb_distance * subgrid$zagreb_dist+
  subcoef$bfe1 * subgrid$krajina6_id1+
  subcoef$bfe2 * subgrid$krajina6_id2+
  subcoef$bfe3 * subgrid$krajina6_id3+
  subcoef$bfe4 * subgrid$krajina6_id4+
  subcoef$bfe5 * subgrid$krajina6_id5+
  subcoef$bfe6 * subgrid$krajina6_id6+
  subcoef$bfe7 * subgrid$krajina6_id7+
  subcoef$quad1 * subgrid$quad1+
  subcoef$quad2 * subgrid$quad2+
  subcoef$Intercept * subgrid$ones



log_existing_zadrugas_yokes<-ggplot()+
  geom_raster(data=subgrid, aes(x =xcoord, y = ycoord, fill=fit1))+
  geom_point(data=subdata, 
             aes_string(x='POINT_X', y='POINT_Y', fill="log_existing_zadrugas_yokes"), 
             size=0.7, 
             colour="gray50", shape=21, na.rm=TRUE)+ #pch=21 allows for borders
  geom_sf(data = boundary, colour = "red", linewidth = 1.5)+
  scale_x_continuous("Longitude") +
  scale_y_continuous("Latitude") +
  theme_bw()+
  theme(axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title=element_text(size=14),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))+
  scale_fill_continuous(high="black",low="white", name="Log Units")+
  ggsn::scalebar(x.min = 18.5, x.max = 19.1,
           #Above you mention how long the scalebar should be
           y.min = 46.2, y.max = 46.4,
           #Above you mention how thick the scalebar should be
           dist = 30, dist_unit = "km",
           transform = T, model = "WGS84",
           location = "topright",
           st.size = 4,
           #Above you have the font size of the numbers below the scalebar
           st.dist =0.4,
           #Above you have the distance between the bar and the text, as a proportion of the y axis.
           height=0.2)
log_existing_zadrugas_yokes<-reposition_legend(log_existing_zadrugas_yokes, 'bottom right')
log_existing_zadrugas_yokes

ggsave(log_existing_zadrugas_yokes,file="./Dropbox/Legacies_Project/Paper/figures/Figure_6b.jpg", 
       height=10, width=20, 
       units = "cm", dpi=300)

